Oscillator decomposition of infant fNIRS data

The functional near-infrared spectroscopy (fNIRS) can detect hemodynamic responses in the brain and the data consist of bivariate time series of oxygenated hemoglobin (oxy-Hb) and deoxygenated hemoglobin (deoxy-Hb) on each channel. In this study, we investigate oscillatory changes in infant fNIRS signals by using the oscillator decompisition method (OSC-DECOMP), which is a statistical method for extracting oscillators from time series data based on Gaussian linear state space models. OSC-DECOMP provides a natural decomposition of fNIRS data into oscillation components in a data-driven manner and does not require the arbitrary selection of band-pass filters. We analyzed 18-ch fNIRS data (3 minutes) acquired from 21 sleeping 3-month-old infants. Five to seven oscillators were extracted on most channels, and their frequency distribution had three peaks in the vicinity of 0.01-0.1 Hz, 1.6-2.4 Hz and 3.6-4.4 Hz. The first peak was considered to reflect hemodynamic changes in response to the brain activity, and the phase difference between oxy-Hb and deoxy-Hb for the associated oscillators was at approximately 230 degrees. The second peak was attributed to cardiac pulse waves and mirroring noise. Although these oscillators have close frequencies, OSC-DECOMP can separate them through estimating their different projection patterns on oxy-Hb and deoxy-Hb. The third peak was regarded as the harmonic of the second peak. By comparing the Akaike Information Criterion (AIC) of two state space models, we determined that the time series of oxy-Hb and deoxy-Hb on each channel originate from common oscillatory activity. We also utilized the result of OSC-DECOMP to investigate the frequency-specific functional connectivity. Whereas the brain oscillator exhibited functional connectivity, the pulse waves and mirroring noise oscillators showed spatially homogeneous and independent changes. OSC-DECOMP is a promising tool for data-driven extraction of oscillation components from biological time series data.


State space models and time series decomposition
OSC-DECOMP is based on state space models, which are statistical models of time series data [1,2]. In state space models, we assume a series of unobserved states x 1 , . . . , x N behind the observed time series data y 1 , . . . , y N , where N is the length of the time series. Specifically, a state space model is a pair of the state model and the observation model where v t and w t are called the system noise and the observation noise, respectively. The state model (1) represents the Markov transition of unobserved states, whereas the observation model (2) describes the generation of observed data depending on the current state. In particular, a state space model of the following form is called a Gaussian linear state space model: where x t ∈ R n , y t ∈ R m , F ∈ R n×n , and H ∈ R m×n . A fundamental task in state space models is to estimate the unobserved states x = {x 1 , · · · , x N } from the observed data y = {y 1 , · · · , y N }. We can regard each state variable x t as a parameter of the distribution p(y t | x t ). The state model (1) is then considered to place a prior on x, and the estimation of x from y in (2) is formulated as a Bayes estimation problem. The posterior distribution p(x t | y 1 , · · · , y t ) and p(x s | y 1 , · · · , y t ) with s < t are called the filtering and smoothing distributions, respectively. For general state space models, several algorithms for filtering and smoothing have been developed, such as the particle filter and the ensemble Kalman filter. In particular, for Gaussian linear state space models, filtering and smoothing are realized using the Kalman filter and Kalman smoother, which require only matrix computations. See [2] for details.
State space models can be used to decompose a given time series into several components. For instance, [3] developed a Bayesian seasonal adjustment method for economic time series. They assumed that a given time series y t consists of a trend component T t , seasonal component S t , and observation noise W t : The trend component T t and seasonal component S t are described by the Gaussian linear state space models. By applying the Kalman smoother, the given time series y t is decomposed into these components: Similarly, [4] proposed a method for decomposing a given time series into trend and AR components.

Kalman filter/smoother algorithm
The algorithm of the Kalman filter and smoother is explained in the following. We assume that the initial distribution for x 1 is given by N(x 1|0 , Σ 1|0 ). For example, if the state model (3) is stationary, then we can naturally take the initial distribution as the stationary distribution, which is Gaussian. Then, all of the conditional distributions p(x s | y 1 , · · · , y t ) are Gaussian. Therefore, they are determined by the mean and the covariance matrix. We define the conditional mean as and the conditional covariance matrix as

Kalman filter
We present the Kalman filter algorithm for the Gaussian linear state space model (3) and (4).
The Kalman filter computes the one-step ahead predictive distribution N(x t|t−1 , Σ t|t−1 ), and the filtering distribution N(x t|t , Σ t|t ) sequentially, in the order From the filtering distribution N(x t−1|t−1 , Σ t−1|t−1 ), the one-step ahead predictive distribution N(x t|t−1 , Σ t|t−1 ) is computed as From the one-step ahead predictive distribution N(x t|t−1 , Σ t|t−1 ) and y t , the filtering distribution N(x t|t , Σ t|t ) is computed as is called the Kalman gain.

Kalman smoother
We present the fixed-interval Kalman smoother algorithm for the Gaussian linear state space model (3) and (4).
The fixed-interval Kalman smoother computes the smoothing distribution N(x s|N , Σ s|N ) sequentially in the order It requires the one-step ahead predictive distribution N(x t|t−1 , Σ t|t−1 ) and the filtering distribution N(x t|t , Σ t|t ), which are obtained using the Kalman filter.
From the smoothing distribution N(x t+1|N , Σ t+1|N ), the smoothing distribution N(x t|N , Σ t|N ) is computed as We note that Σ t|N does not depend on y. In particular, it is time-reversible: Σ t|N = Σ N −t+1|N .

Hessian computation and confidence interval
Here, we present an extension of the Kalman filter algorithm for computing the gradient and Hessian of the log-likelihood of the Gaussian linear state space models. It is a generalization of the algorithm in [5] to multi-dimensional observations. We also explain the method for constructing confidence intervals by using the computed Hessian.
Consider a Gaussian linear state space model where x t ∈ R n , y t ∈ R m , F ∈ R n×n , and H ∈ R m×n . In the following, we omit the argument θ for simplicity. The one-step ahead predictive distribution of y t is where ε t and R t are the one-step ahead prediction error and its covariance defined by Therefore, the log-likelihood is given by

Gradient computation
By differentiating (6), the gradient of the log-likelihood is where, from (5), From the Kalman filter formula, the required derivatives are computed sequentially as and where K t = Σ t|t−1 H R −1 t is the Kalman gain.

Hessian computation
By differentiating (7) and using the Hessian of the log-likelihood is where, from (9), From (10)-(14), the required derivatives are computed sequentially as and